Molecular Dynamics#

Fig. 28 MD vs MC: Both sample microstates. The former follows the natural motion (dynamics) the later samples from Botlzman distribution following rules designed to make sampling efficient.#
Timescales and Lengthscales#
Classical Molecular Dynamics can access a hiearrchy of time-scales from pico seconds to microseconds.
It is also possible to go beyond the time scale of brute force MD byb emplying clever enhanced sampling techniques.
Fig. 29 Different time-scales underlying different leng-scales/motions in molecules#

Is MD just Newton’s laws applied on big systems?#
Not quite: Noble prize in Chemistry 2013
Classical molecular dynamics (MD) is a powerful computational technique for studying complex molecular systems.
Applications span wide range including proteins, polymers, inorganic and organic materials.
Alos molecular dynamics simulation is being used in a complimentary way to the analysis of experimental data coming from NMR, IR, UV spectroscopies and elastic-scattering techniques, such as small angle scattering or diffraction.
Integrating equations of motion#
The simplest integrating scheme for ODEs is the Euler’s method.
Given the \(n\)-dimensional vectors from the ODE standard form we are writing down equation in finite difference form.
Much better integrators are known under the names of Runge Kutta, 2nd, 4th, 6th … order
import numpy as np
import matplotlib.pyplot as plt
def euler(y, f, t, h):
"""Euler integrator: Returns new y at t+h.
"""
return y + h * f(t, y)
def rk2(y, f, t, h):
"""Runge-Kutta RK2 midpoint"""
k1 = f(t, y)
k2 = f(t + 0.5*h, y + 0.5*h*k1)
return y + h*k2
def rk4(y, f, t, h):
"""Runge-Kutta RK4"""
k1 = f(t, y)
k2 = f(t + 0.5*h, y + 0.5*h*k1)
k3 = f(t + 0.5*h, y + 0.5*h*k2)
k4 = f(t + h, y + h*k3)
return y + h/6 * (k1 + 2*k2 + 2*k3 + k4)
Integrate Harmonic Oscillator using Euler, RK2, RK4
Show code cell source
def f(t, y):
''' Define a simple harmonic potential'''
return np.array([y[1], -y[0]])
y = np.array([1., 1.])
pos, vel = [], []
t = 0
h = 0.1
for i in range(1000):
y = euler(y, f, t, h) # Change integration method venv, euler, rk2, rk4
t+=h
pos.append(y[0])
vel.append(y[1])
fig, ax = plt.subplots(ncols=3, nrows=1, figsize=(12,3))
pos, vel =np.array(pos), np.array(vel)
ax[0].plot(pos, vel)
ax[1].plot(pos)
ax[1].plot(vel)
ax[2].plot(0.5*pos**2 + 0.5*vel**2)
[<matplotlib.lines.Line2D at 0x7f42bc13f340>]
Verlet algortihm#
Taylor expansion of position \(\vec{r}(t)\) after timestep \(\Delta t\) we obtain forward and backward Euler schems
In 1967 Loup Verlet introduced a new algorithm into molecular dynamics simulations which preserves energy is accurate and efficient.
Summing the two taylor expansion above we get a updating scheme which is an order of mangnitude more accurate
Velocity is not needed to update the positions. But we still need them to set the temperature.
Terms of order \(O(\Delta t^3)\) cancel in position giving position an accuracy of order \(O(\Delta t^4)\)
To update the position we need positions in the past at two different time points! This is is not very efficient.
Velocity Verlet updating scheme#
A better updating scheme has been proposed known as Velocity-Verlet (VV) which stores positions, velocities and accelerations at the same time. Time stepping backward expansion \(r(t-\Delta t + \Delta t)\) and summing with the forward Tayloer expansions we get Velocity Verlet updating scheme:
Substituting forces \(a=\frac{F}{m}\) instead of acelration we get
Velocity Verlet Algorithm
1. Evaluate the initial force from the current position:
2. Update the position:
3. Partially update the velocity:
4. Evaluate the force at the new position:
5. Complete the velocity update:
def velv(y, f, t, h):
"""Velocity Verlet for solving differential equations.
A little inefficient since same force (first and last) is evaluated twice!
"""
# 1. Evluate force
F = f(t, y)
# 2, Velocity partial update
y[1] += 0.5*h * F[1]
# 3. Full step position
y[0] += h*y[1]
# 4. Force re-eval
F = f(t+h, y)
# 5. Full step velocity
y[1] += 0.5*h * F[1]
return y
Show code cell source
def f(t, y):
''' Define a simple harmonic potential'''
return np.array([y[1], -y[0]])
y = np.array([1., 1.])
pos, vel = [], []
t = 0
h = 0.1
for i in range(1000):
y = velv(y, f, t, h) # Change integration method venv, euler, rk2, rk4
t+=h
pos.append(y[0])
vel.append(y[1])
fig, ax = plt.subplots(ncols=3, nrows=1, figsize=(12,3))
pos, vel =np.array(pos), np.array(vel)
ax[0].plot(pos, vel)
ax[1].plot(pos)
ax[1].plot(vel)
ax[2].plot(0.5*pos**2 + 0.5*vel**2)
[<matplotlib.lines.Line2D at 0x7f42b7ec11c0>]
Ensemble averages#
Ergodic hypothesis
\(\tau_{eq}\) time it teks for system to “settle into” equilibrium and \(t\) time of simulation.
Time averages are equal to ensemble averages when \(t\lim \infty\). But in real life we need to determine how long to sample.
Ergodic assumption is: during the coruse of MD we visit all the microstates (or small but representative subset) that go into the ensemble average!
Energy conservation in simulations
Kinetic, potential and total energies are flcutuation quantities in the Molecualr dynamics simulations.
Energy in (NVE) or its average \(\langle E \rangle\) (in NVT, NPT, etc) must remain constant!
Temperature control in simulations
According to equipariting result of equilibrium statistical mechanics in the NVT ensmeble
Pressure control in simulations
Molecular Dynamics of Classical Harmonic Oscillator (NVE)#
Define the Verlet integrator to loop over steps and record data#
Show code cell source
import numpy as np
def run_md_nve_1d(x, v, dt, t_max, en_force):
"""
Minimalistic 1D Molecular Dynamics simulation (NVE ensemble)
using Velocity Verlet integration.
Simulates a particle moving in a 1D potential without thermal noise or friction
(i.e., energy-conserving dynamics).
Parameters
----------
x : float
Initial position.
v : float
Initial velocity.
dt : float
Time step for integration.
t_max : float
Total simulation time.
en_force : callable
Function that takes position `x` as input and returns a tuple (potential energy, force).
Returns
-------
pos : ndarray
Array of particle positions over time.
vel : ndarray
Array of particle velocities over time.
KE : ndarray
Array of kinetic energies over time.
PE : ndarray
Array of potential energies over time.
Example
-------
>>> def harmonic_force(x):
>>> k = 1.0
>>> return 0.5 * k * x**2, -k * x
>>> pos, vel, KE, PE = run_md_nve_1d(1.0, 0.0, 0.01, 10.0, harmonic_force)
"""
times, pos, vel, KE, PE = [], [], [], [], []
# Initialize force and potential energy
pe, F = en_force(x)
t = 0.0
for step in range(int(t_max / dt)):
# Velocity Verlet Integration
# Half-step velocity update
v += 0.5 * F * dt
# Full-step position update
x += v * dt
# Update force at new position
pe, F = en_force(x)
# Half-step velocity update
v += 0.5 * F * dt
# Save results
times.append(t)
pos.append(x)
vel.append(v)
KE.append(0.5 * v * v)
PE.append(pe)
# Advance time
t += dt
return np.array(pos), np.array(vel), np.array(KE), np.array(PE)
Apply to a potential of harmonic oscillator#
Show code cell source
#----parameters of simulation----
k = 3
x0 = 1
v0 = 0
dt = 0.01 * 2*np.pi/np.sqrt(k) #A good timestep determined by using oscillator frequency
t_max = 1000
### Define Potential Energy function
def ho_en_force(x, k=k):
'''Force field of harmonic oscillator:
returns potential energy and force'''
return k*x**2, -k*x
### Run simulation
pos, vel, KE, PE = run_md_nve_1d(x0, v0, dt, t_max, ho_en_force)
### Visualize
fig, ax =plt.subplots(ncols=2)
ax[0].hist(pos);
ax[1].hist(vel, color='orange');
ax[0].set_ylabel('P(x)')
ax[1].set_ylabel('P(v)')
fig.tight_layout()
Langevin equation#
A particle of mass \(m\) moves under the force derived from a potential energy \(U(x)\). The motion is purely deterministic.
Challenge: What should we do when we have only a one or few particles, and cannot explicitly simulate the vast surrounding environment in order to assign temperature?
Solution: We model the surrounding medium (e.g., a solvent) as an implicit thermal bath that interacts with the particle.
The particle exchanges energy with the bath, maintaining thermal equilibrium at a fixed temperature \(T\).
This motivates Langevin dynamics, where the effects of the solvent are captured by friction (dissipation) and random thermal kicks (fluctuations), without simulating solvent molecules explicitly.
Langevin Equation
Overdamped Limit of Langevin Dynamics (\(m \ddot{x} = 0\))
The friction \(\lambda\) and thermal noise \(\eta(t)\) are clearly connected becasue the faster the particle movies (more noise) the more it also dissipates energy.
the connection is known by the name of Fluctuation-Dissipation Theorem:
Fluctuation-Dissipation Theorem (FDT)
The environment “forgets” what happened almost immediately after a collision very short memory hence noise terms are uncorrelated (independent). Just like what we had in brownian motion.
The FDT ensures that the strength of random thermal kicks is precisely tuned to the amount of viscous damping, so that the system reaches and maintains thermal equilibrium at temperature \(T\)
FDT guarantees that thermal fluctuations (random noise) and viscous dissipation are balanced, preventing runaway motion and stabilizing equilibrium.
FDT connects diffusion (random spreading) and viscosity (resistance to motion), both fundamentally controlled by temperature.
Derivation of FDT
Let us focus on the overdamped Langevin equation for simplicity:
Assume a flat potential first (i.e., \(U(x) = 0\)) to understand the role of friction and noise alone. (You can later generalize to nonzero \(U(x)\), but this captures the key idea.)
Thus, integrating over time:
Let’s compute the mean-square displacement:
Expanding:
Now — if we assume \(\eta(t)\) is delta-correlated (white noise):
where \(C\) is some constant to be determined. Then:
This shows diffusive behavior with an effective diffusion constant:
because for normal diffusion we expect:
Connect to equilibrium (Einstein relation)
From statistical mechanics, the equilibrium distribution for \(x\) must obey the equipartition theorem:
If we define the equilibrium probability distribution \(P(x)\) for the free particle:
(no restoring force). But the velocity (or here the “rate of change”) still carries thermal energy:
In equilibrium, for a particle moving slowly through a viscous medium, we know:
The diffusion coefficient \(D\) satisfies the Einstein relation:
Thus, matching:
Multiplying both sides by \(2\lambda^2\), we find:
This is the fluctuation-dissipation theorem for Langevin dynamics:
Fluctuations (noise strength) are linked to dissipation (friction \(lambda\)) and temperature \(T\)).
Summary Intuition
Without friction, noise would make energy grow indefinitely → no equilibrium.
Without noise, friction would freeze the system.
FDT guarantees just the right balance to maintain thermal equilibrium at temperature \(T\).
Molecualr Dynamics of Harmonic oscillator (NVT)#
Our goal is to Numerically solve the Langevin equation:
The target is to sample the canonical distribution:
There are many ways of simulating langevin we will use BAOAB scheme.
The main idea is to Split the Langevin evolution operator into three pieces:
B: Free drift in position (\(x\) evolution under \(v\)).
A: Velocity change under conservative force \(F = -\nabla V(x)\).
O: Stochastic Ornstein-Uhlenbeck (OU) process (friction + random kicks in \(v\)).
Thus, BAOAB updates position and velocity via a B → A → O → A → B sequence at each timestep \(\Delta t\).
Feature |
Explanation |
|---|---|
Stability |
Handles both large and small friction \(\gamma\) robustly |
Accuracy |
Samples correct Boltzmann distribution up to \(\mathcal{O}(\Delta t^2)\) errors |
Efficiency |
Large stable timesteps allowed compared to naive Euler schemes |
Universality |
Reduces to velocity Verlet (Hamiltonian MD) if \(\gamma=0\) and to overdamped Langevin if \(\gamma \to \infty\) |
BAOAB scheme
B: Drift half-step in position
\[ x \to x + \frac{\Delta t}{2} v \]A: Kick velocity half-step using force
\[ v \to v + \frac{\Delta t}{2m} F(x) \]O: Apply stochastic friction and noise to velocity
\[ v \to e^{-\gamma \Delta t} v + \sqrt{k_B T (1 - e^{-2\gamma \Delta t})/m} \, R, \]where (R) is a standard normal random number.
A: Kick velocity half-step again
\[ v \to v + \frac{\Delta t}{2m} F(x) \]B: Drift half-step in position
\[ x \to x + \frac{\Delta t}{2} v \]
Show code cell source
def langevin_md_1d(x, v, dt, kBT, gamma, t_max, en_force):
'''Langevin dynamics applied to 1D potentials
Using integration scheme known as BA-O-AB.
INPUT: Any 1D function with its parameters
'''
times, pos, vel, KE, PE = [], [], [], [], []
t = 0
for step in range(int(t_max/dt)):
#B-step
pe, F = en_force(x)
v += F*dt/2
#A-step
x += v*dt/2
#O-step
v = v*np.exp(-gamma*dt) + np.sqrt(1-np.exp(-2*gamma*dt)) * np.sqrt(kBT) * np.random.normal()
#A-step
x += v*dt/2
#B-step
pe, F = en_force(x)
v += F*dt/2
### Save output
times.append(t), pos.append(x), vel.append(v), KE.append(0.5*v*v), PE.append(pe)
return np.array(times), np.array(pos), np.array(vel), np.array(KE), np.array(PE)
Run lagevin dynamics of 1d harmonic oscillator#
Show code cell source
import numpy as np
import matplotlib.pyplot as plt
# Initial conditions
x0 = 0.0
v0 = 0.5
# Input parameters
kBT = 0.25
gamma = 0.01
dt = 0.01
t_max = 10000
freq = 10
k = (2 * np.pi * freq)**2 # Spring constant (harmonic oscillator)
### Define Potential: Energy and Force
def ho_en_force(x, k=k):
energy = 0.5 * k * x**2
force = -k * x
return energy, force
### Langevin dynamics (assuming you have this function correctly defined)
# Should return arrays: times, pos, vel, KE, PE
times, pos, vel, KE, PE = langevin_md_1d(x0, v0, dt, kBT, gamma, t_max, ho_en_force)
### Plotting
fig, ax = plt.subplots(nrows=1, ncols=4, figsize=(10, 4))
bins = 50
# Theoretical distributions
def gaussian_x(x, k, kBT):
return np.exp(-k*x**2/(2*kBT)) / np.sqrt(2*np.pi*kBT/k)
def gaussian_v(v, kBT):
return np.exp(-v**2/(2*kBT)) / np.sqrt(2*np.pi*kBT)
# Plot histograms
ax[0].hist(pos, bins=bins, density=True, alpha=0.6, color='skyblue')
ax[1].hist(vel, bins=bins, density=True, alpha=0.6, color='salmon')
# Plot theoretical curves
x_grid = np.linspace(min(pos), max(pos), 300)
v_grid = np.linspace(min(vel), max(vel), 300)
ax[0].plot(x_grid, gaussian_x(x_grid, k, kBT), 'k-', lw=2, label='Theory')
ax[0].set_xlabel('Position x')
ax[0].set_ylabel('P(x)')
ax[0].legend()
ax[1].plot(v_grid, gaussian_v(v_grid, kBT), 'k-', lw=2, label='Theory')
ax[1].set_xlabel('Velocity v')
ax[1].set_ylabel('P(v)')
ax[1].legend()
ax[2].plot(pos[-1000:], vel[-1000:], alpha=0.3)
ax[2].set_xlabel('Position x')
ax[2].set_ylabel('Velocity v')
ax[2].set_title('Phase Space')
E_tot=KE+PE
ax[3].plot(E_tot)
ax[3].set_xlabel('Time')
ax[3].set_ylabel('Energy')
ax[3].set_title('Total Energy')
fig.tight_layout()
plt.show()
Show code cell source
import numpy as np
import plotly.graph_objects as go
# Compute total energy
KE = 0.5 * vel**2
E = KE + PE
# Subsample for performance (optional)
stride = 2
x = pos[::stride]
v = vel[::stride]
E = (KE + PE)[::stride]
# Create the 3D scatter plot
fig = go.Figure(data=go.Scatter3d(
x=x, y=v, z=E,
mode='markers',
marker=dict(
size=2,
color=E, # Color by total energy
colorscale='Viridis',
colorbar=dict(title='E'),
opacity=0.8
)
))
fig.update_layout(
title='Langevin Trajectory in (x, v, E) Space',
scene=dict(
xaxis_title='Position x',
yaxis_title='Velocity v',
zaxis_title='Total Energy E'
),
width=800,
height=600
)
fig.show()